#remettre ligne de code general identification (trop long a calcul pour les test de knit) + remettre pour cluster 13 deux premiere ligne list..
#cluster 3 DP small, 8,10,16
#11 refaire text marqueur
#refaire cluster 13 avec annot one note pour petit iddentif des pop
#DP blast put cell cycle umap ? dans le fichier myc_pten_paper --> report mycpten..cccaintegration_modif
#revoir titre partie en fonction du papier
#refaire docker en verifiant que tout les librairies necessaire sont dowloand dedans
#verif tout comm sont bein en anglais
#CD8 vs CD4 remettre les calcul de fin marker pour le final et pas le load fichier qui est la poru gagner du temps
#loading final object obtain with Experiment_preprocessing
SAMPLE1 <- "181031"
SAMPLE2 <- "190211"
if(! file.exists(paste0(OUTPUT_PATH, "T-Seurat-merged_clean-subset",".Robj"))){
print("You should start with Experiment_preprocessing.Rmd or dowload our final object 'T-Seurat-merged_clean-subset.Robj' ")
do <- FALSE
}else{
print ("You are starting analysis of our final Seurat object")
load(paste0(OUTPUT_PATH, "T-Seurat-merged_clean-subset",".Robj"))
do <- TRUE
}
[1] “You are starting analysis of our final Seurat object”
Idents(T.Seurat) <- "HTO"
T.Seurat.thymus <- subset(T.Seurat, idents = c("Myc- PTEN- thymus","MYC- thymus","PTEN- thymus","WT thymus"))
a <- t(margin.table(table(T.Seurat.thymus@meta.data$HTO,T.Seurat.thymus@meta.data$integrated_snn_res.1.8),2))
T.Seurat.spleen <- subset(T.Seurat, idents = c("Myc- PTEN- spleen","MYC- spleen","PTEN- spleen","WT spleen"))
b <- t(margin.table(table(T.Seurat.spleen@meta.data$HTO,T.Seurat.spleen@meta.data$integrated_snn_res.1.8),2))
c <- t((a/(a+b)*100))
c
##
## [,1]
## 0 8.204812
## 1 16.918967
## 2 99.730942
## 3 99.634369
## 4 2.474794
## 5 1.868132
## 6 77.661431
## 7 95.885510
## 8 99.564270
## 9 3.139013
## 10 38.990826
## 11 3.240741
## 12 1.960784
## 13 15.763547
## 14 82.543641
## 15 9.595960
## 16 44.074074
## 17 37.327189
## 18 98.913043
## 19 0.000000
## 20 94.771242
## 21 99.264706
## 22 80.882353
## 23 100.000000
#Thymic populations
thymus.clusters <- rownames(as.data.frame(c[which(c[,1]>25),]))
Idents(T.Seurat.thymus) <- "integrated_snn_res.1.8"
T.Seurat.thymus <- subset(T.Seurat.thymus, idents = thymus.clusters)
DimPlot(T.Seurat.thymus)
#Spleen populations
as.data.frame(100-c[which(c[,1]<75),])
## 100 - c[which(c[, 1] < 75), ]
## 0 91.79519
## 1 83.08103
## 4 97.52521
## 5 98.13187
## 9 96.86099
## 10 61.00917
## 11 96.75926
## 12 98.03922
## 13 84.23645
## 15 90.40404
## 16 55.92593
## 17 62.67281
## 19 100.00000
spleen.clusters <- rownames(as.data.frame(c[which(c[,1]<75),]))
Idents(T.Seurat.spleen) <- "integrated_snn_res.1.8"
T.Seurat.spleen <- subset(T.Seurat.spleen, idents = spleen.clusters)
DimPlot(T.Seurat.spleen)
Based on several markers we are adjusting our clusterging, similar cluster are then annotate as one.
Thymic clusters :
#Look at differentiation markers on thymic clusters
DotPlot(T.Seurat.thymus, features = c("percent.mito","Bmf","Trp53inp1","Tox2","Cd5","Cd69","Cd27","Rag1","Rag2","Cd4","Cd8a","Cd8b1","Uchl3","Mki67","Cdk1","Ptcra","Il2ra","Cd34")) + scale_colour_gradient2(low = "steelblue", mid = "white", high = "red") + theme_dark(base_size = 14) + coord_flip()
#Regroup thymic clusters
#Similar cluster are annotate as one
Idents(T.Seurat.thymus) <- "integrated_snn_res.1.8"
T.Seurat.thymus@meta.data$manualclusters = "nothing"
T.Seurat.thymus@meta.data[WhichCells(T.Seurat.thymus, slot = "integrated_snn_res.1.8", idents = "22"),]$manualclusters = "22"
T.Seurat.thymus@meta.data[WhichCells(T.Seurat.thymus, slot = "integrated_snn_res.1.8", idents = "21"),]$manualclusters = "21"
T.Seurat.thymus@meta.data[WhichCells(T.Seurat.thymus, slot = "integrated_snn_res.1.8", idents = c("20","18","14")),]$manualclusters = "20,18,14"
T.Seurat.thymus@meta.data[WhichCells(T.Seurat.thymus, slot = "integrated_snn_res.1.8", idents = c("8","2","3","23")),]$manualclusters = "8,2,3,23"
T.Seurat.thymus@meta.data[WhichCells(T.Seurat.thymus, slot = "integrated_snn_res.1.8", idents = "6"),]$manualclusters = "6"
T.Seurat.thymus@meta.data[WhichCells(T.Seurat.thymus, slot = "integrated_snn_res.1.8", idents = "7"),]$manualclusters = "7"
T.Seurat.thymus@meta.data[WhichCells(T.Seurat.thymus, slot = "integrated_snn_res.1.8", idents = c("10","16")),]$manualclusters = "10,16"
T.Seurat.thymus@meta.data[WhichCells(T.Seurat.thymus, slot = "integrated_snn_res.1.8", idents = "17"),]$manualclusters = "17"
Splenic clusters :
#Look at differentiation markers on splenic clusters
DotPlot(T.Seurat.spleen, features = c("Aes","Anxa1","Ifng","Itgal","Foxp3","S1pr1","Sell","Ccr7","Trdc","Tcrg-C4","Tcrg-C2","Tcrg-C1","Trbc2","Trbc1","Trac","Cd8b1","Cd4")) + scale_colour_gradient2(low = "steelblue", mid = "white", high = "red") + theme_dark(base_size = 14) + coord_flip()
#Regroup splenic clusters
#Similar cluster are annotate as one
Idents(T.Seurat.spleen) <- "integrated_snn_res.1.8"
T.Seurat.spleen@meta.data$manualclusters = "nothing"
T.Seurat.spleen@meta.data[WhichCells(T.Seurat.spleen, slot = "integrated_snn_res.1.8", idents = "4"),]$manualclusters = "4"
T.Seurat.spleen@meta.data[WhichCells(T.Seurat.spleen, slot = "integrated_snn_res.1.8", idents = "1"),]$manualclusters = "1"
T.Seurat.spleen@meta.data[WhichCells(T.Seurat.spleen, slot = "integrated_snn_res.1.8", idents = "9"),]$manualclusters = "9"
T.Seurat.spleen@meta.data[WhichCells(T.Seurat.spleen, slot = "integrated_snn_res.1.8", idents = "12"),]$manualclusters = "12"
T.Seurat.spleen@meta.data[WhichCells(T.Seurat.spleen, slot = "integrated_snn_res.1.8", idents = "11"),]$manualclusters = "11"
T.Seurat.spleen@meta.data[WhichCells(T.Seurat.spleen, slot = "integrated_snn_res.1.8", idents = "19"),]$manualclusters = "19"
T.Seurat.spleen@meta.data[WhichCells(T.Seurat.spleen, slot = "integrated_snn_res.1.8", idents = "13"),]$manualclusters = "13"
T.Seurat.spleen@meta.data[WhichCells(T.Seurat.spleen, slot = "integrated_snn_res.1.8", idents = "17"),]$manualclusters = "17"
T.Seurat.spleen@meta.data[WhichCells(T.Seurat.spleen, slot = "integrated_snn_res.1.8", idents = c("10","16")),]$manualclusters = "10,16"
T.Seurat.spleen@meta.data[WhichCells(T.Seurat.spleen, slot = "integrated_snn_res.1.8", idents = "15"),]$manualclusters = "15"
T.Seurat.spleen@meta.data[WhichCells(T.Seurat.spleen, slot = "integrated_snn_res.1.8", idents = c("0","5")),]$manualclusters = "0,5"
Setting thymic cluster order
Idents(T.Seurat.thymus) <- "manualclusters"
T.Seurat.thymus@active.ident <- factor(T.Seurat.thymus@active.ident,levels=c("22","21","20,18,14","8,2,3,23","7","6","10,16","17"))
levels(T.Seurat.thymus)
## [1] "22" "21" "20,18,14" "8,2,3,23" "7" "6"
## [7] "10,16" "17"
DotPlot(T.Seurat.thymus, dot.scale = 8,features = c("percent.mito","Bmf","Trp53inp1","Tox2","Cd5","Cd69","Cd27","Rag1","Rag2","Cd4","Cd8a","Cd8b1","Mki67","Cdk1","Ptcra","Il2ra","Cd34")) + scale_colour_gradient2(low = "steelblue", mid = "white", high = "red") + coord_flip()+ theme(
panel.background = element_rect(fill = "grey65",
size = 0.5, linetype = "solid"),
panel.grid.major = element_line(size = 0.5, linetype = 'solid',
colour = "grey55"))
Setting splenic cluster order
Idents(T.Seurat.spleen) <- "manualclusters"
T.Seurat.spleen@active.ident <- factor(T.Seurat.spleen@active.ident,levels=c("10,16","1","4","0,5","12","15","9","17","13","11","19"))
DotPlot(T.Seurat.spleen, dot.scale = 8, features= c("Foxp3","Aes","Anxa1","Gzma","Ccl5","Cxcr3","Sell","S1pr1","Ccr7","Trdc","Tcrg-C4","Tcrg-C2","Tcrg-C1","Trbc2","Trbc1","Trac","Cd8b1","Cd4")) + scale_colour_gradient2(low = "steelblue", mid = "white", high = "red") + theme(
panel.background = element_rect(fill = "grey65",
size = 0.5, linetype = "solid"),
panel.grid.major = element_line(size = 0.5, linetype = 'solid',
colour = "grey55")) + coord_flip()
Idents(T.Seurat) <- "integrated_snn_res.1.8"
DefaultAssay(T.Seurat) <- "RNA"
#markerspleen <- FindAllMarkers(T.Seurat)
Idents(T.Seurat) <- "integrated_snn_res.1.8"
FeaturePlot(T.Seurat, features = c("adt_CD4","Sell","Cd8b1"),ncol = 3,cols = c("grey", "light blue","cyan3","cyan4","dodgerblue3","blue","mediumslateblue","purple","orchid3","red","brown","black"))
FeaturePlot(T.Seurat, features = c("Cd4","Cd8b1"),blend = T,
blend.threshold= 0.4,cols = c("grey", "light blue","cyan3","cyan4","dodgerblue3","blue","mediumslateblue","purple","orchid3","red","brown","black"))
Cluster 0 seems to contain naive SP4 cells (CD4+, Cd8- and CD62L-)
FeaturePlot(T.Seurat, features = c("adt_CD4","Cd4","Cd8b1","Sell"),ncol = 4,cols = c("grey", "light blue","cyan3","cyan4","dodgerblue3","blue","mediumslateblue","purple","orchid3","red","brown","black"))
FeaturePlot(T.Seurat, features = c("Cd4","Cd8b1"),blend = T,
blend.threshold= 0.4,cols = c("grey", "light blue","cyan3","cyan4","dodgerblue3","blue","mediumslateblue","purple","orchid3","red","brown","black"))
Cluster 1 seems to contain CD8 naive T cells (CD4-, Cd8+ and CD62L-)
FeaturePlot(T.Seurat, features = c("adt_CD4","Cd4","Cd8b1","Sell"),ncol = 4,cols = c("grey", "light blue","cyan3","cyan4","dodgerblue3","blue","mediumslateblue","purple","orchid3","red","brown","black"))
FeaturePlot(T.Seurat, features = c("Cd4","Cd8b1"),blend = T,
blend.threshold= 0.4,cols = c("grey", "light blue","cyan3","cyan4","dodgerblue3","blue","mediumslateblue","purple","orchid3","red","brown","black"))
Cluster 2 seems to contain DP cells (CD4+, CD8+, CD62L-)
FeaturePlot(T.Seurat, features = c("adt_CD4","Cd4","Cd8b1","Sell"),ncol = 4,cols = c("grey", "light blue","cyan3","cyan4","dodgerblue3","blue","mediumslateblue","purple","orchid3","red","brown","black"))
FeaturePlot(T.Seurat, features = c("Cd4","Cd8b1"),blend = T,
blend.threshold= 0.4,cols = c("grey", "light blue","cyan3","cyan4","dodgerblue3","blue","mediumslateblue","purple","orchid3","red","brown","black"))
Cluster 4 seems to contain CD8 naive T cells (CD4-, Cd8+, and CD62L-)
FeaturePlot(T.Seurat, features = c("adt_CD4","Sell","Cd8b1"),ncol = 3,cols = c("grey", "light blue","cyan3","cyan4","dodgerblue3","blue","mediumslateblue","purple","orchid3","red","brown","black"))
FeaturePlot(T.Seurat, features = c("Cd4","Cd8b1"),blend = T,
blend.threshold= 0.4,cols = c("grey", "light blue","cyan3","cyan4","dodgerblue3","blue","mediumslateblue","purple","orchid3","red","brown","black"))
Cluster 5 seems to contain naive SP4 cells (CD4+, Cd8- and CD62L-)
FeaturePlot(T.Seurat, features = c("percent.mito"),cols = c("grey", "light blue","cyan3","cyan4","dodgerblue3","blue","mediumslateblue","purple","orchid3","red","brown","black"), pt.size = 1)
Cluster 6 seems to contain dying cells
FeaturePlot(T.Seurat, features = c("Cd69","Cd5","Satb1"),cols = c("grey", "light blue","cyan3","cyan4","dodgerblue3","blue","mediumslateblue","purple","orchid3","red","brown","black")) #marker of TCR activation
Seems to contain cells going from DP to SP between thymus and spleen
FeaturePlot(T.Seurat, features = c("Cd4","adt_CD4"))
VlnPlot(T.Seurat.spleen, features = c("Sell","Ccr7","Il7r","Cxcr3"),ncol=2) # Cxcr3 high on effector
Seems to contain CD4 eff (Cd4+, Sell-, Ccr7-, Cd127-)
FeaturePlot(T.Seurat, features = c("Cd8b1","Sell","Ccl5","adt_CD4"),cols = c("grey", "light blue","cyan3","cyan4","dodgerblue3","blue","mediumslateblue","purple","orchid3","red","brown","black"))
FeaturePlot(T.Seurat, features = c("Trac","Trbc1","Trdc","Tcrg-C1"),cols = c("grey", "light blue","cyan3","cyan4","dodgerblue3","blue","mediumslateblue","purple","orchid3","red","brown","black"))
FeaturePlot(T.Seurat, features = c("Ly6c2","Nkg7"),cols = c("grey", "light blue","cyan3","cyan4","dodgerblue3","blue","mediumslateblue","purple","orchid3","red","brown","black"))
VlnPlot(T.Seurat.spleen, features = c("Sell","Ccr7","Il7r"))
Cluster 11 Seems to contain CD8 mem cells & seem to contain Tgd cells & Ly6c2+
FeaturePlot(T.Seurat, features = c("adt_CD4","Sell","Cd8b1"),cols = c("grey", "light blue","cyan3","cyan4","dodgerblue3","blue","mediumslateblue","purple","orchid3","red","brown","black"))
FeatureScatter(object = T.Seurat, feature1 = "adt_CD4", feature2 = "CD8", cells = colnames(subset(T.Seurat, idents = "0")), slot = "data")
Cluster 12 seems to contain naive SP4 cells (CD4+, Cd8- and CD62L -)
plot1 <- DotPlot(cluster13, dot.scale = 10, features= c("Foxp3","Aes","Anxa1","Gzma","Ccl5","Cxcr3","Sell","S1pr1","Ccr7","Trdc","Tcrg-C4","Tcrg-C2","Tcrg-C1","Trbc2","Trbc1","Trac","Cd8b1","Cd4")) + scale_colour_gradient2(low = "steelblue", mid = "white", high = "red") + theme(
panel.background = element_rect(fill = "grey65",
size = 0.5, linetype = "solid"),
panel.grid.major = element_line(size = 0.5, linetype = 'solid',
colour = "grey55")) + coord_flip()
plot2 <- DotPlot(cluster13, dot.scale = 10, features= unique(top_genes_feature_plot$gene)) + scale_colour_gradient2(low = "steelblue", mid = "white", high = "red") + theme(
panel.background = element_rect(fill = "grey65",
size = 0.5, linetype = "solid"),
panel.grid.major = element_line(size = 0.5, linetype = 'solid',
colour = "grey55")) + coord_flip()
grid.arrange(plot1, plot2, ncol=2)
#table(cluster13@meta.data$integrated_snn_res.1)
This cluster is heterogeneous
FeaturePlot(T.Seurat, features = c("adt_CD4","Cd4","Cd8b1"),cols = c("grey", "light blue","cyan3","cyan4","dodgerblue3","blue","mediumslateblue","purple","orchid3","red","brown","black"))
VlnPlot(T.Seurat, features = c("Pcna","Cdk1","Mki67")) #teichmann
Cluster 14 seems to be Dp blast
FeaturePlot(T.Seurat, features = c("adt_CD4","Cd4","Foxp3","Il2ra"),cols = c("grey", "light blue","cyan3","cyan4","dodgerblue3","blue","mediumslateblue","purple","orchid3","red","brown","black"))
Cluster 15 seems to contain Treg (CD4+, Cd8-, Cd25+, Foxp3+)
FeaturePlot(T.Seurat, features = c("adt_CD4","Cd8b1","Trdc","Tcrg-C1"),cols = c("grey", "light blue","cyan3","cyan4","dodgerblue3","blue","mediumslateblue","purple","orchid3","red","brown","black"))
Cluster 17 seem to contain Tgd DN (CD4-, Cd8-, TCR D+, TCR G+)
FeaturePlot(T.Seurat, features = c("adt_CD4","Cd4","Cd8b1"),cols = c("grey", "light blue","cyan3","cyan4","dodgerblue3","blue","mediumslateblue","purple","orchid3","red","brown","black"))
VlnPlot(T.Seurat, features = c("Pcna","Cdk1","Mki67")) #teichmann
Dp blast
FeaturePlot(T.Seurat, features = c("Cd8b1","Sell","Ccl5","Gzma","Klrc1","Ifng","Zeb2","Ly6c2"),ncol = 3,cols = c("grey", "light blue","cyan3","cyan4","dodgerblue3","blue","mediumslateblue","purple","orchid3","red","brown","black"))
# Zeb2 terminally differentiated CTL
VlnPlot(T.Seurat.spleen, features = c("Sell","Ccr7","Il7r"))
Cluster 19 seemsto contain CD8 CTL (CD4-, Cd8+, Ccl5+, Sell-, ccr7-,cd127-) & Ly6c2+
FeaturePlot(T.Seurat, features = c("adt_CD4","Cd4","Cd8b1"),cols = c("grey", "light blue","cyan3","cyan4","dodgerblue3","blue","mediumslateblue","purple","orchid3","red","brown","black"))
VlnPlot(T.Seurat, features = c("Pcna","Cdk1","Mki67")) #teichmann
Seems to be Dp blast
FeaturePlot(T.Seurat, features = c("Cd8b1","adt_CD4"))
VlnPlot(T.Seurat, features = c("Cd8b1","adt_CD4"))
Cluster 21 seem to contain ISP cells (CD4-, Cd8+)
FeaturePlot(T.Seurat, features = c("adt_CD4","Cd4","Cd8b1","Il2ra"),cols = c("grey", "light blue","cyan3","cyan4","dodgerblue3","blue","mediumslateblue","purple","orchid3","red","brown","black"))
Cluster 22 seems to contain DN T cells (CD4-, Cd8-, Cd25+)
VlnPlot(T.Seurat, features = c("Cd8b1","adt_CD4"))
VlnPlot(T.Seurat, features = c("Rag1")) #on DP cells mostly (immgen), DP quicient or DP blast (Teichman)
Seems to be Dp small (Cd4+,Cd8+,Rag1+)
#check genotype proportion in each spleen clusters
df <- as.data.frame(as.data.frame.matrix(t(prop.table(table(T.Seurat.spleen@meta.data$HTO,T.Seurat.spleen@meta.data$manualclusters),1)*100)))
df$cluster = rownames(df)
df2<-as.data.frame(t(cbind(rep(60,11),rep(0,11),df)[,1:6]))
rownames(df2[1:2,]) <- c("60","0")
#order data frame
df2 <- df2[c("10,16","19","11","13","17","9","15","12","0,5","4","1")]
radarchart(df2, cglcol="grey", cglty=1 ,cglwd=0.8, vlcex=0.8, pcol=c("#FF99FF","coral1","cyan3","chartreuse3") , plwd=3, plty=1 ,caxislabels=paste(seq(from = 0,to = 60,by = 15),"%"), axislabcol = "grey40", axistype = 0)
#check genotype proportion in each thymic cluster
df <- as.data.frame(as.data.frame.matrix(t(prop.table(table(T.Seurat.thymus@meta.data$HTO,T.Seurat.thymus@meta.data$manualclusters),1)*100)))
df$cluster = rownames(df)
df2<-as.data.frame(t(cbind(rep(65,8),rep(0,8),df)[,1:6]))
rownames(df2[1:2,]) <- c("65","0")
#order
df2 <- df2[c("22","17","10,16","6","7","8,2,3,23","20,18,14","21")]
radarchart(df2, cglcol="grey", cglty=1,caxislabels=paste(seq(0,70,17.5),"%"), axistype = 0,axislabcol="grey40", cglwd=0.8, vlcex=0.8, pcol=c("#FF99FF","coral1","cyan3","chartreuse3") , plwd=3, plty=1 )
### Proportion bar plot
#thymic proportions
datathym <- data.frame(prop.table(t(prop.table(table(T.Seurat.thymus@meta.data$HTO,T.Seurat.thymus@meta.data$manualclusters),1)),1)*100)
#order cluster level
datathym$Var1 <- factor(datathym$Var1, levels = c("22","21","20,18,14","8,2,3,23","7","6","10,16","17"))
cellnumber <- data.frame(colSums(table(T.Seurat.thymus@meta.data$HTO,T.Seurat.thymus@meta.data$manualclusters)) )
cellnumber$cluster <- rownames(cellnumber)
row_order <- c("22","21","20,18,14","8,2,3,23","6","7","17","10,16")
cellnumber <- cellnumber[row_order,]
Cluster <- datathym$Var1
HTO <- datathym$Var2
Percentage <- datathym$Freq
text <- cellnumber$colSums.table.T.Seurat.thymus.meta.data.HTO..T.Seurat.thymus.meta.data.manualclusters..
cols <- c("Myc- PTEN- thymus" = "#FF99FF", "MYC- thymus" = "coral1", "PTEN- thymus" = "cyan3", "WT thymus" = "chartreuse3")
ggplot(datathym, aes(fill=HTO, y=Percentage, x=Cluster)) +
geom_bar(position="stack", stat="identity", width=0.7)+
xlab("Thymic Clusters")+ylab("Percentage")+ theme(legend.position="bottom")+ scale_fill_manual(values =cols, labels=c("Myc Pten","Myc","Pten","WT"))+theme_light()+geom_hline(yintercept=c(25,50,75), linetype="dashed", color = "grey60")+ annotate("text", x = c(1,2,3,4,5,6,7,8), y=103, label = c(paste0(text)))
#splenic proportions
Idents(T.Seurat.spleen) <- "manualclusters"
T.Seurat.spleenbar <- subset(T.Seurat.spleen, idents = c("0,5","9","1","11","19")) #to only keep cluster needed for the plot
data3 <- data.frame(prop.table(t(prop.table(table(T.Seurat.spleenbar@meta.data$HTO,T.Seurat.spleenbar@meta.data$manualclusters),1)),1)*100)
#order cluster level
data3$Var1 <- factor(data3$Var1, levels = c("0,5","9","1","11","19"))
cellnumber <- data.frame(colSums(table(T.Seurat.spleen@meta.data$HTO,T.Seurat.spleen@meta.data$manualclusters)) )
cellnumber$cluster <- rownames(cellnumber)
row_order <-c("0,5","9","1","11","19")
cellnumber <- cellnumber[row_order,]
Cluster <- data3$Var1
HTO <- data3$Var2
Percentage <- data3$Freq
text <- cellnumber$colSums.table.T.Seurat.spleen.meta.data.HTO..T.Seurat.spleen.meta.data.manualclusters..
# Stacked bar plot
cols <- c("Myc- PTEN- spleen" = "#FF99FF", "MYC- spleen" = "coral1", "PTEN- spleen" = "cyan3", "WT spleen" = "chartreuse3")
ggplot(data3, aes(fill=HTO, y=Percentage, x=Cluster)) +
geom_bar(position="stack", stat="identity", width=0.7)+
xlab("Splenic Clusters")+ylab("Percentage")+ theme(legend.position="bottom")+scale_fill_manual(values =cols, labels=c("Myc Pten","Myc","Pten","WT")) +theme_light()+geom_hline(yintercept=c(25,50,75), linetype="dashed", color = "grey60")+ annotate("text", x = c(1,2,3,4,5), y=103, label = c(text))
Idents(T.Seurat) <- "integrated_snn_res.1.8"
#DGE between our two CD8 naive clusters
dgecd8_data <- FindMarkers(T.Seurat, ident.1 = "1", ident.2 = "4",logfc.threshold=0)
# bar plot
#add abs value to table
dgecd8_data$abs <- abs(dgecd8_data$avg_logFC)
dgecd8_data$genename <- rownames(dgecd8_data)
# Select markers for plotting on a Heatmap
markers.use=subset(dgecd8_data, p_val_adj<1e-50 & abs>0.20)
dfcd8markers <-markers.use[order(markers.use$avg_logFC),]
dfcd8markers$genename <- factor(dfcd8markers$genename, levels = dfcd8markers$genename[order(dfcd8markers$avg_logFC)])
dfcd8markers$logpval <- log10(dfcd8markers$p_val_adj)
ggplot(dfcd8markers, aes(x = dfcd8markers$genename, y = dfcd8markers$avg_logFC, fill = logpval)) + # Fill column
geom_bar(stat = "identity", width = .6) + # draw the bars
ylim(-1.2,1.2)+
labs(title="Try dGE") +
theme_tufte() + # Tufte theme from ggfortify
theme(plot.title = element_text(hjust = .5),axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
axis.ticks = element_blank()) +
scale_fill_gradient2(low='red', mid='orange', high='blue',midpoint = -120, breaks=c(-52,-120,-200),labels=c("-50","-120","-200"))+coord_flip() # Flip axes
#### CD4
#define proportion of genotype in CD4 naive clusters
T.Seurat.spleen@meta.data$manualHTO = "nothing"
Idents(T.Seurat.spleen) <- "MULTI_ID"
T.Seurat.spleen@meta.data[WhichCells(T.Seurat.spleen, slot = "MULTI_ID", idents = c("Spleen-ctrl","Spleen-P")),]$manualHTO = "Spleen-ctrl&P"
T.Seurat.spleen@meta.data[WhichCells(T.Seurat.spleen, slot = "MULTI_ID", idents = c("Spleen-M","Spleen-MP")),]$manualHTO = "Spleen-M&MP"
prop.table(t(prop.table(table(T.Seurat.spleen@meta.data$integrated_snn_res.1.8,T.Seurat.spleen@meta.data$manualHTO),1)),2)*100
##
## 0 1 2 3 4 5 6 7 8
## Spleen-ctrl&P 17.473118 20.257235 94.642857 28.667413
## Spleen-M&MP 82.526882 79.742765 5.357143 71.332587
##
## 9 10 11 12 13 14
## Spleen-ctrl&P 75.000000 25.187970 33.732057 71.000000 60.818713
## Spleen-M&MP 25.000000 74.812030 66.267943 29.000000 39.181287
##
## 15 16 17 18 19 20 21 22 23
## Spleen-ctrl&P 43.016760 36.423841 17.647059 35.256410
## Spleen-M&MP 56.983240 63.576159 82.352941 64.743590
#in our CD4 naive clusters (0,5 and 12) the one with most cells from M&MP is 0 and the one with most cells from ctrl&P is 12
Idents(T.Seurat) <- "integrated_snn_res.1.8"
#DGE between 0 and 12
dgecd4_data <- FindMarkers(T.Seurat, ident.1 = "0", ident.2 = "12",logfc.threshold=0)
# bar plot
#add abs value to table
dgecd4_data$abs <- abs(dgecd4_data$avg_logFC)
dgecd4_data$genename <- rownames(dgecd4_data)
# Select markers for plotting on a Heatmap
markers.use=subset(dgecd4_data,p_val_adj<1e-10 & abs>0.2)
dfcd4markers <-markers.use[order(markers.use$avg_logFC),]
dfcd4markers$genename <- factor(dfcd4markers$genename, levels = dfcd4markers$genename[order(dfcd4markers$avg_logFC)])
dfcd4markers$logpval <- log10(dfcd4markers$p_val_adj)
ggplot(dfcd4markers, aes(x = dfcd4markers$genename, y = dfcd4markers$avg_logFC, fill = logpval)) + # Fill column
geom_bar(stat = "identity", width = .6) + # draw the bars
ylim(-1.2,1.2)+
labs(title="Try dGE CD4") +
theme_tufte() + # Tufte theme from ggfortify
theme(plot.title = element_text(hjust = .5),axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
axis.ticks = element_blank()) +
scale_fill_gradient2(low='red', mid='orange', high='blue',midpoint = -40, breaks=c(-12,-40,-57),labels=c("-10","-40","-60"))+ coord_flip() # Flip axes
Idents(T.Seurat) <- "integrated_snn_res.1.8"
# get all gene name express in our cells as background
background <- T.Seurat@assays$RNA@meta.features
backgroundrow <- rownames(background)
#DGE between our two CD8 naive clusters
#dgecd8_data <- FindMarkers(T.Seurat, ident.1 = "1", ident.2 = "4",logfc.threshold=0)
#upreg_cd8 <- subset(dgecd8_data, avg_logFC>0)
#downreg_cd8 <- subset(dgecd8_data, avg_logFC<0.2 & p_val_adj<1e-50)
#genecomprow <- rownames(downreg_cd8)
#### TO SUPPRESS FOR FINAL
#write.table(genecomprow,file="/home/nozaism/Workspace/01_These/01_Project/Myc_Pten_Paper/cluster1v4.txt",sep="\t",quote=F)
genecomprow <- read.table("/home/nozaism/Workspace/01_These/01_Project/Myc_Pten_Paper/cluster1v4.txt", sep = "\t")
genecomprow$x = as.character(genecomprow$x)
genecomprow <- genecomprow[,1]
####
CPenrich <- enrichGO(gene= genecomprow, OrgDb = 'org.Mm.eg.db', ont="BP",keyType = "SYMBOL",universe = backgroundrow) # org.Mm.eg.db genome mouse
head (CPenrich)
## ID Description GeneRatio
## GO:0002181 GO:0002181 cytoplasmic translation 16/51
## GO:0042254 GO:0042254 ribosome biogenesis 21/51
## GO:0022613 GO:0022613 ribonucleoprotein complex biogenesis 22/51
## GO:0042274 GO:0042274 ribosomal small subunit biogenesis 13/51
## GO:0042255 GO:0042255 ribosome assembly 13/51
## GO:0000028 GO:0000028 ribosomal small subunit assembly 8/51
## BgRatio pvalue p.adjust qvalue
## GO:0002181 79/13400 2.587342e-24 1.495484e-21 1.342694e-21
## GO:0042254 263/13400 4.286807e-23 1.238887e-20 1.112314e-20
## GO:0022613 397/13400 9.509083e-21 1.832083e-18 1.644904e-18
## GO:0042274 61/13400 3.841622e-20 5.551144e-18 4.983999e-18
## GO:0042255 66/13400 1.180651e-19 1.364832e-17 1.225391e-17
## GO:0000028 17/13400 5.863805e-16 5.648799e-14 5.071677e-14
## geneID
## GO:0002181 Rplp0/Rpl18a/Rps28/Rpl17/Rpl39/Rps2/Rps23/Rpl10a/Rpl15/Rplp1/Rpl8/Rpl18/Rpl36/Rpl24/Rps29/Rpl26
## GO:0042254 Rps19/Rplp0/Rps5/Rps28/Rps24/Rps7/Rps6/Rps2/Rpl23a/Rps16/Rpl3/Rpsa/Rps8/Rpl10a/Rpl12/Rps14/Rpl14/Rps10/Rps15/Rpl24/Rpl26
## GO:0022613 Rps19/Rplp0/Rps5/Rps28/Rps24/Rps7/Rps6/Rps2/Rpl23a/Rps16/Rpl3/Rpsa/Rps8/Rps23/Rpl10a/Rpl12/Rps14/Rpl14/Rps10/Rps15/Rpl24/Rpl26
## GO:0042274 Rps19/Rps5/Rps28/Rps24/Rps7/Rps6/Rps2/Rps16/Rpsa/Rps8/Rps14/Rps10/Rps15
## GO:0042255 Rps19/Rplp0/Rps5/Rps28/Rps2/Rpl23a/Rpl3/Rpsa/Rpl12/Rps14/Rps10/Rps15/Rpl24
## GO:0000028 Rps19/Rps5/Rps28/Rps2/Rpsa/Rps14/Rps10/Rps15
## Count
## GO:0002181 16
## GO:0042254 21
## GO:0022613 22
## GO:0042274 13
## GO:0042255 13
## GO:0000028 8
dotplot(CPenrich, showCategory=15,color = "p.adjust",x="count") #+ coord_flip()+theme(axis.text.x = element_text(angle = 90, hjust = 1))
emapplot(CPenrich, showCategory = 50)
#DGE between our two most control and Myd del CD4 naive clusters
#dgecd4_data <- FindMarkers(T.Seurat, ident.1 = "0", ident.2 = "12",logfc.threshold=0)
#downreg_cd4 <- subset(dgecd4_data, avg_logFC<0& p_val_adj<1e-20)
#genecomprow <- rownames(downreg_cd4)
#### TO SUPPRESS FOR FINAL
#write.table(genecomprow,file="/home/nozaism/Workspace/01_These/01_Project/Myc_Pten_Paper/cluster0vs12.txt",sep="\t",quote=F)
genecomprow <- read.table("/home/nozaism/Workspace/01_These/01_Project/Myc_Pten_Paper/cluster1v4.txt", sep = "\t")
genecomprow$x = as.character(genecomprow$x)
genecomprow <- genecomprow[,1]
####
CPenrich <- enrichGO(gene= genecomprow, OrgDb = 'org.Mm.eg.db', ont="BP",keyType = "SYMBOL",universe = backgroundrow) # org.Mm.eg.db genome mouse
dotplot(CPenrich, showCategory=15,color = "p.adjust",x="count")
Deciphering what are our eYFP negative cells on effector and memory CD8 T cells
Heatmap :
# On cluster 11 - cd8 memory
C11Ctrl <-rownames(T.Seurat.spleen@meta.data[T.Seurat.spleen@meta.data$MULTI_ID == "Spleen-ctrl" & T.Seurat.spleen@meta.data$integrated_snn_res.1.8 == "11",])
mC11Ctrl <- mean(T.Seurat.spleen@assays$RNA@data["eYFP",C11Ctrl])
C11Pt <- rownames(T.Seurat.spleen@meta.data[T.Seurat.spleen@meta.data$MULTI_ID == "Spleen-P" & T.Seurat.spleen@meta.data$integrated_snn_res.1.8 == "11",])
mC11Pt <- mean(T.Seurat.spleen@assays$RNA@data["eYFP",C11Pt])
C11MP <- rownames(T.Seurat.spleen@meta.data[T.Seurat.spleen@meta.data$MULTI_ID == "Spleen-MP" & T.Seurat.spleen@meta.data$integrated_snn_res.1.8 == "11",])
mC11MP <- mean(T.Seurat.spleen@assays$RNA@data["eYFP",C11MP])
C11M <- rownames(T.Seurat.spleen@meta.data[T.Seurat.spleen@meta.data$MULTI_ID == "Spleen-M" & T.Seurat.spleen@meta.data$integrated_snn_res.1.8 == "11",])
mC11M <- mean(T.Seurat.spleen@assays$RNA@data["eYFP",C11M])
#on cluster 19 - cd8 eff term
C19Ctrl <-rownames(T.Seurat.spleen@meta.data[T.Seurat.spleen@meta.data$MULTI_ID == "Spleen-ctrl" & T.Seurat.spleen@meta.data$integrated_snn_res.1.8 == "19",])
mC19Ctrl <- mean(T.Seurat.spleen@assays$RNA@data["eYFP",C19Ctrl])
C19Pt <- rownames(T.Seurat.spleen@meta.data[T.Seurat.spleen@meta.data$MULTI_ID == "Spleen-P" & T.Seurat.spleen@meta.data$integrated_snn_res.1.8 == "19",])
mC19Pt <- mean(T.Seurat.spleen@assays$RNA@data["eYFP",C19Pt])
C19MP <- rownames(T.Seurat.spleen@meta.data[T.Seurat.spleen@meta.data$MULTI_ID == "Spleen-MP" & T.Seurat.spleen@meta.data$integrated_snn_res.1.8 == "19",])
mC19MP <- mean(T.Seurat.spleen@assays$RNA@data["eYFP",C19MP])
C19M <- rownames(T.Seurat.spleen@meta.data[T.Seurat.spleen@meta.data$MULTI_ID == "Spleen-M" & T.Seurat.spleen@meta.data$integrated_snn_res.1.8 == "19",])
mC19M <- mean(T.Seurat.spleen@assays$RNA@data["eYFP",C19M])
#on cluster 1
C1Ctrl <-rownames(T.Seurat.spleen@meta.data[T.Seurat.spleen@meta.data$MULTI_ID == "Spleen-ctrl" & T.Seurat.spleen@meta.data$integrated_snn_res.1.8 == "1",])
mC1Ctrl <- mean(T.Seurat.spleen@assays$RNA@data["eYFP",C1Ctrl])
C1Pt <- rownames(T.Seurat.spleen@meta.data[T.Seurat.spleen@meta.data$MULTI_ID == "Spleen-P" & T.Seurat.spleen@meta.data$integrated_snn_res.1.8 == "1",])
mC1Pt <- mean(T.Seurat.spleen@assays$RNA@data["eYFP",C1Pt])
C1MP <- rownames(T.Seurat.spleen@meta.data[T.Seurat.spleen@meta.data$MULTI_ID == "Spleen-MP" & T.Seurat.spleen@meta.data$integrated_snn_res.1.8 == "1",])
mC1MP <- mean(T.Seurat.spleen@assays$RNA@data["eYFP",C1MP])
C1M <- rownames(T.Seurat.spleen@meta.data[T.Seurat.spleen@meta.data$MULTI_ID == "Spleen-M" & T.Seurat.spleen@meta.data$integrated_snn_res.1.8 == "1",])
mC1M <- mean(T.Seurat.spleen@assays$RNA@data["eYFP",C1M])
l1 <- c("1","1","1","1","11","11","11","11","19","19","19","19")
l2 <- c("Ctrl","Pten","Myc","MycPten","Ctrl","Pten","Myc","MycPten","Ctrl","Pten","Myc","MycPten")
l3 <- c(mC1Ctrl,mC1Pt,mC1M,mC1MP,mC11Ctrl,mC11Pt,mC11M,mC11MP,mC19Ctrl,mC19Pt,mC19M,mC19MP)
tableauYFP <- data.frame(Cluster = l1, Genotypes =l2)
tableauYFP <- cbind(tableauYFP,Values = l3)
tableauYFP$Cluster <- factor(tableauYFP$Cluster,levels = c("1","11","19"))
tableauYFP$Genotypes <- factor(tableauYFP$Genotypes,levels = c("MycPten","Myc","Pten","Ctrl"))
ggplot(tableauYFP, aes(x = Cluster, Genotypes)) +
geom_tile(aes(fill = Values)) +
scale_fill_gradient2( mid='yellow', high='red',limits=c(0,max(tableauYFP$Values)))+ theme_classic(base_size=20)
eYFP neagtive are coming from Myc and Myc Pten del mice
#DOT plot eyFP and TGD on CD8 effector and memory
Idents(T.Seurat.spleen) <- "integrated_snn_res.1.8"
CD8sub <- subset(T.Seurat.spleen, idents = c("11","19","13"))
Idents(CD8sub) <- "MULTI_ID"
CD8sub@active.ident <- factor(CD8sub@active.ident,levels=c("Spleen-M","Spleen-MP","Spleen-ctrl","Spleen-P"))
levels(CD8sub)
## [1] "Spleen-M" "Spleen-MP" "Spleen-ctrl" "Spleen-P"
DotPlot(CD8sub, dot.scale = 8,features = c("Tcrg-C1","Trdc","Trbc2","Trac","eYFP") ) + scale_colour_gradient2(low = "steelblue", mid = "white", high = "red")+ ggtitle("CD8 memory and effector clusters")
eyFP negative are coming from Myc and Myc Pten del mice and are Tgd cells.
sessionInfo()
## R version 3.5.3 (2019-03-11)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.5 LTS
##
## Matrix products: default
## BLAS/LAPACK: /usr/lib/libopenblasp-r0.2.18.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=C
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] org.Mm.eg.db_3.7.0 AnnotationDbi_1.44.0 IRanges_2.16.0
## [4] S4Vectors_0.20.1 Biobase_2.42.0 BiocGenerics_0.28.0
## [7] knitr_1.23 RColorBrewer_1.1-2 magrittr_1.5
## [10] dplyr_0.8.1 ggthemes_4.2.0 clusterProfiler_3.10.1
## [13] enrichplot_1.2.0 fmsb_0.6.3 gridExtra_2.3
## [16] kableExtra_1.1.0 plotly_4.9.0 ggplot2_3.1.1
## [19] Seurat_3.0.1
##
## loaded via a namespace (and not attached):
## [1] fastmatch_1.1-0 plyr_1.8.4 igraph_1.2.4.1
## [4] lazyeval_0.2.2 splines_3.5.3 BiocParallel_1.16.6
## [7] listenv_0.7.0 urltools_1.7.3 digest_0.6.19
## [10] htmltools_0.3.6 GOSemSim_2.8.0 viridis_0.5.1
## [13] GO.db_3.7.0 gdata_2.18.0 memoise_1.1.0
## [16] cluster_2.0.9 ROCR_1.0-7 globals_0.12.4
## [19] readr_1.3.1 graphlayouts_0.6.0 R.utils_2.8.0
## [22] prettyunits_1.0.2 colorspace_1.4-1 blob_1.1.1
## [25] rvest_0.3.4 ggrepel_0.8.1 xfun_0.7
## [28] crayon_1.3.4 jsonlite_1.6 survival_2.44-1.1
## [31] zoo_1.8-5 ape_5.3 glue_1.3.1
## [34] polyclip_1.10-0 gtable_0.3.0 webshot_0.5.1
## [37] UpSetR_1.4.0 future.apply_1.2.0 scales_1.0.0
## [40] DOSE_3.8.2 DBI_1.0.0 bibtex_0.4.2
## [43] Rcpp_1.0.1 metap_1.1 viridisLite_0.3.0
## [46] progress_1.2.2 gridGraphics_0.5-0 reticulate_1.12
## [49] bit_1.1-14 europepmc_0.3 rsvd_1.0.0
## [52] SDMTools_1.1-221.1 tsne_0.1-3 htmlwidgets_1.3
## [55] httr_1.4.0 fgsea_1.8.0 gplots_3.0.1.1
## [58] ica_1.0-2 pkgconfig_2.0.2 R.methodsS3_1.7.1
## [61] farver_2.0.3 labeling_0.3 ggplotify_0.0.5
## [64] tidyselect_0.2.5 rlang_0.3.4 reshape2_1.4.3
## [67] munsell_0.5.0 tools_3.5.3 RSQLite_2.1.1
## [70] ggridges_0.5.1 evaluate_0.13 stringr_1.4.0
## [73] yaml_2.2.0 npsurv_0.4-0 bit64_0.9-7
## [76] fitdistrplus_1.0-14 tidygraph_1.1.2 caTools_1.17.1.2
## [79] purrr_0.3.2 RANN_2.6.1 ggraph_2.0.2
## [82] pbapply_1.4-0 future_1.13.0 nlme_3.1-140
## [85] R.oo_1.22.0 DO.db_2.9 xml2_1.2.0
## [88] compiler_3.5.3 rstudioapi_0.10 png_0.1-7
## [91] lsei_1.2-0 tibble_2.1.1 tweenr_1.0.1
## [94] stringi_1.4.3 lattice_0.20-38 Matrix_1.2-17
## [97] pillar_1.4.0 BiocManager_1.30.4 Rdpack_0.11-0
## [100] triebeard_0.3.0 lmtest_0.9-37 data.table_1.12.2
## [103] cowplot_0.9.4 bitops_1.0-6 irlba_2.3.3
## [106] gbRd_0.4-11 qvalue_2.14.1 R6_2.4.0
## [109] KernSmooth_2.23-15 codetools_0.2-16 MASS_7.3-51.4
## [112] gtools_3.8.1 assertthat_0.2.1 withr_2.1.2
## [115] sctransform_0.2.0 hms_0.4.2 grid_3.5.3
## [118] tidyr_0.8.3 rmarkdown_1.12 rvcheck_0.1.8
## [121] Rtsne_0.15 ggforce_0.3.1 base64enc_0.1-3